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The integral method can be used to model accurately flows down an inclined plane. Such 
a method consists in projecting the full 3D equations on a lower dimensional representa- 
tion. The vertical velocity profiles have their functional form fixed, based from the exact 
solution of homogeneous steady flows. This projection is achieved by integration of the 
momentum equation over the flow depth - Saint- Venant approach. Here we generalize 
the viscous case to two non-newtonian constitutive relations: a Prandtl-like turbulent 
closure and a local granular rheology. We discuss one application in each case: the forma- 
tion of anti-dunes in viscous streams, the transverse velocity profile in turbulent channels 
and the Kapitza instability in dense granular flows. They demonstrate the usefulness of 
this approach to get a model qualitatively correct, quantitatively reasonable and in which 
the dynamical mechanisms at work can be easily identified. 



1. Introduction 

The accurate hydrodynamical description of free surface flows, although initiated 
during the XIX th , still remains a modern problem of physics. The difficulties to be 
overcome are of different sorts. In some situations, like a viscous flow on a slope, 
the governing equations are both well posed and easy to integrate numerically. How- 
ever, a numerical experiment, although giving access to any measurement, by con- 
trast to actual experiments, does not provide in itself a physical understanding of the 
system. Most often, theories aiming to enlighten the physical mechanisms at work 
have been based on low dimensional model equations - like the lubrication approxi- 
mation for viscous flows ([Landau fc Lifschitz 1959)) . In other situations, the numeri- 
cal inegration of the governing equations is not accessible, due to the range of rele- 
vant length-scales and time-scales involved. For instance, the simulation of turbulent 
flows is currently limited to slightly less than four decades in space while the simula- 
tion of a river section would require at least 10 4 times more mesh-points. Even for 
viscous flows, the presence of a moving contact line implies a divergence of the vis- 
cous stress regularised at the molecular scale. Five to six decades of length-scales being 
involved, the simulation of a simple millimetre scale moving droplet is beyond usual 
computer facilities. One trick is to introduce a sub-grid model, like for Large Eddy 
Simulations of turbulence. Another solution consists again in deriving a low dimen- 
sional model. In the turbulent case, one can use the shallow water equations with- 
out (|Landau fc Lifschitz 1959) IWhitham 1974)) or with a damping term (|Julien 20021) . 
In the viscous case, the lubrication approximation can be extended to moving contact 
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lines by introducing the slip length flSnoeijer et al. 2006} |Snoeijer et al. 20071 ) . Finally, 
there exist situations, like dense granular flows, for which the constitutive relations have 
not been derived from first principles and are still under debate (jGDR MIDI 2004)) . 
In the particular case of granular matter, the depth averaged dynamical equation 
hereafter referred as the Saint- Venant equation (|St-Ven ant 1843j ISt-Venant 1850[) - has 
been widely used to bypass the problem. Within this framework, instead of prescrib- 
ing a stress/strain relation, one simply models the friction on the static bed over which 
grains are flowin g (| Savage fe Hutter 19891 |Savage fe Hutter 199T] IDouady et al. 1999] 
Pouliquen 1999bl|Pouliquen fe Forterre 20021|Mangeney et al. 2003) . 



In any flow bounded by a free surface h(x,y,t), the idea is to simplify the full set of 
equations governing the evolution of the velocity u(x, y, z,t) (3D equations) into that 
for the depth averaged velocity U(x,y,t) (referred to as (2+l)D equations). To perform 
this reduction under controlled approximations, one usually assumes that the flow is 
thin, i.e. that the longitudinal variations are negligible in comparison to transverse 
ones: 0(d x ) <C 0(d z ). This restricts the validity of the approach to small interface 
slopes, i.e. \d x h\ <C 1. To go beyond this restriction, methods in the spirit of the 
lubrication approximation have been developped to recover quantitative results for larger 
slopes. For instance, a long- wavelength expansion has been performed by |Snoeijer 2 006 
for wedge-like geometries, where the restriction is on the curvature of the interface, i.e. 
\hd xx h\ <C 1. It yields to an equation that is well adapted to treat contact line problems 
dSnoeijer et al. 20071IDelon et al. 2"007] . 

However, by construction, the lubrication approximation (even extended to larger 
slopes) is lacking important dynamical mechanisms, namely inertial effects and longi- 
tudinal diffusion of momentum. It is therefore unable to give much insight into in- 
stability and pattern formation mechanisms. This is why a lot of numerical models, 
starting from a long-wave expansion, have been derived since the pioneering work of 
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the fluid into large-amplitude strongly non-linear localized structures (jPumir et al. 19 83). 
Assuming that the height of the waves remains small in comparison to their wave- 
length, i.e. \h — (h)\ -C A, the interface modulations are expected to be governed by 
a Bcnney-like equation (Bcnncy 1966), where the temporal evolution of h is expressed 
as a function of various algebraic powers of h and its gradients. The most recent models 
( |Ruyer-Quil fc M anncville 1998, 20001 120041) accurately and economically predict linear 
and non-linear properties of film flows up to relatively high Reynolds numbers. 

Such developments do not exist for real turbulent or granular flows. In this paper, 
we show how one can, in a general way, take into account all terms of Navier-Stokes 
equations, while keeping the Saint- Venant projection. This goal is achieved by looking 
for approximate solutions (called below 'test functions') whose form is inferred from that 
of the exact steady and homogeneous solutions of the original equations, and with which 
we require that these equations are fulfilled on average only, i.e. integrated over the flow 
depth. Of course, we expect this approach to work well for situations which are not too 
far from the steady and homogeneous case. 

We shall explain the precise development of this integral method in the next section 
for the case of the well studied viscous flow. We then turn to more original calculations 
with the treatment of turbulent and granular flows. For each of these three cases, we 
illustrate our approach with a concrete example of application: the formation of dunes 
in viscous streams, the transverse velocity profile in turbulent channels and the Kapitza 
instability in dense granular flows. Although we are greatly interested in the physical 
insights of these different situations, we shall here focus on the technical discussion of 
the benefits of this approach. 



Integral method for flows down an incline: viscous, turbulent and granular cases 3 




Figure 1. Schematic drawing of the flow down an inclined of slope tan#. x is along the bottom, 
z perpendicular to it. y is in the third direction. Z(x) is the equation of the bed and H(x) that 
of the free surface. We note R = H — Z the local fluid depth, q is the local flow rate - see 
equation (|2,5[l - so that U = q/R is the mean velocity at position x. 

2. Viscous flow equations 

2.1. Governing equations 

We consider the generic situation of a viscous fluid flowing down a plane inclined at 
an angle 9 with respect to the horizontal. The fluid motion is referred to a Cartesian 
coordinate system in which the x-axis coincides with the inclined bottom. The z-axis is 
perpendicular to this plane and y is in the plane, transverse to the slope. A schematic 
drawing of the flow with all notations is shown in figure [TJ Combined with the incom- 
pressibility equation, the Navier-Stokes equation provides the theoretical framework that 
describes the motion of the fluid: 

d jUj = 0, 

d t Ui + (ujdj)ui = -dip + g { + vl\u u (2.1) 

where u — (u, v, w) is the velocity field, p is the pressure divided by the fluid density 
p, g the gravity vector and v the kinematic viscosity. Note that in the entire paper, we 
assume the density to be constant and omit this 1/p factor for simplicity. However, we 
shall still use the word 'pressure' for p - same for momentum and stress. The boundary 
conditions we wish to impose are the following: 

(i) all velocity components vanish at the bottom: Ui(z = Z) = 0, 

(ii) the free surface follows the fluid motion, 

d t H + ud x H + vd y H-w = 0, (2.2) 

all velocity components u, v and w being evaluated at z = H , 
(hi) the stress normal to the surface vanishes. 

2.2. Exact homogeneous steady solutions 

The strategy of the integral method proposed here is to solve exactly the case of ho- 
mogeneous steady solutions and to use it as a test function to derive the Saint- Venant 
equations. In the viscous case, equation (|2.1[) admits a trivial solution corresponding to 
a steady constant-thickness film. Assuming that all derivatives with respect to t, x and 
y vanish, there is only a streamwise velocity u and the equations reduce to: 



= v d zz u + g sin (9, 
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= -d z p- gcosO, (2.3) 

which, for a film of thickness Hq (in this homogeneous case, Zq = and thus Rq = Hq), 
yields to 

g sin 9 

uo(z) = 2 ^ z(2H a - z), 

Po{z) = Patm + 5(^0 - z) COS 0, (2.4) 

where p a t m is the reference pressure. 

2.3. Integral method 

The next step is to reduce the full equations to a set of two equations governing the 
evolution of the local depth of the flow R = H — Z . For the sake of simplicity we shall 
keep to 1 + 1 dimensions, but the following derivation can be easily generalized to (2+l)D. 
The local flow rate q is defined as: 

i-H 

udz. (2.5) 

z 

From q, we define the mean velocity U — qj R. The first of these equations is simply 
the (integral) mass conservation: dtR + d x q — 0. The second one is the projection of 
the Navier-Stokes equation integrated over the height of the flow. It gives an integral 
momentum conservation equation (Saint- Venant equation): 

/ [dtu + ud x u + wd z u] dz = I [g sin 9 — d x p + vA 
J z J z 



+ vAu] dz. (2.6) 



In order to compute these integrals, the vertical profiles are assumed to remain close to 
those obtained in the homogeneous steady case. We thus take a parabolic velocity test 
profile parametrized by q, R and Z (or H) 

» = 3l(i^Vl-^V (2.7) 



R\ R J \ 2R 

Similarly, the pressure profile remains hydrostatic: 

P = Patm+ g(H - z) cos 6». (2.8) 

We shall turn to the limit of this assumption in the conclusion. 

Because of the incompressibility relation, we can compute the z-component of the 
velocity: 

d x udz. (2.9) 

z 

An analysis of the order of magnitude of this velocity leads, as expected, to a ratio w/u ~ 
d x Z . In the following we shall then focus on the streamwise component velocity only, 
as we are interested in situations with a weak bed slope correction to the homogeneous 
Z(x) = case. Using incompressibility as well as the free surface and bottom boundary 
conditions, one can write the inertial terms as the derivative of the momentum flux: 

[d t u + ud x u + wd z u] dz = d t udz + d x u 2 dz = d t q + d x a — I . (2.10) 
z Jz Jz V R J 

Interestingly, the value of a depends only on the test function and is thus determined in 
the homogeneous steady case: a — 6/5 for the parabolic profile (|2.7p . 
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The right hand side of the Saint- Venant equation contains the different forces inte- 
grated vertically, namely: 

(i) the projection of gravity along the streamwise direction, gRsin9, 

(ii) the pressure gradient, — gRcos 9d x H, 

(iii) the vertical diffusion of momentum, 



r-H 



d z {yd z u) dz = ~vd z u\ z= ^ 



o 9 



(iv) the streamwise diffusion of momentum 
d x (vd x u) dz = vd xx q — 3v 



We finally obtain the following system of equations: 



Z + R = H, 
d t R + d x q = 



(2.11) 



(2.12) 



(2.13) 



d t q + d x 



R 



— gRcos 8 (tan (9 
d x (|) d x 



d x H) - ( 1 + -Rd xx H ) + vd xx q, 



3v 



H 



R? 



(d x zy 



(2.14) 



In the long wavelength approximation, at low Reynolds number, one can neglect, in 
equation (|2.14|1 . its left-hand side as well as the terms coming from (I2.12p . One then 
recovers the standard lubrication equations. In the right-hand side of equation (|2.12p . 
one recognizes a term corresponding to the viscous diffusion of q. The three other dis- 
sipative terms originate from the geometry of the flow and their interpretation is less 
straightforward . 

The above equation set is very close in spirit to the derivation of Shkado v 19671 Other 
versions of this modeling have been proposed for specific problems (Ru yer-Quil fc Mam icvillc 1998, 
12000) l2004|) . They all keep the same functional form, but with different numerical co- 
efficients, adjusted to match experimental data. By contrast, our equations are less 
quantitative but directly come from the original Navier-Stokes equations and the dif- 
ferent terms have therefore not been modified in an ad hoc way for some particular 
purpose. All dynamical mechanisms are thus present and can be backtracked to get 
their physical meaning. The interest of the terms extending the lubrication approxima- 
tion for the problem of Kapitza instability and roll waves has been discussed at length 
by | Ruyer-Quil fc Manneville 1998") We refer the interested reader to the review section 
of that article. In the following we discuss another example: the formation of sand 
anti-dunes in a viscous stream. 

2.4. Steady flow over a bumpy bottom 

Let us first apply our approach to the simple case of a steady flow over a non-uniform 
bed Z(x). From conservation of mass, we infer that the flow rate is uniform along the 
stream and is equal to: 

1 gH$ sin 9 



q = H U = 



(2.15) 



3 v 

The flow may be usefully characterized by two independent non-dimensional numbers 
among the following three: the slope tan 9, the Froude number 

JJ 2 



Fr z 



gH cos 9 



(2.16) 




0.01 0.1 1 10 Fr 100 



Figure 2. (a) Amplitude (solid line) and phase (dashed line) of the free surface with respect 
to the bed deformations as a function of the Froude number, for tan# = 0.1 and k = tv/2. The 
free surface profiles have been plotted for a low Froude number (b) and for a high one (c). 



and the Reynolds number 

U H 5 Fr 2 



Re 



v 2 tan 8 

In the limit of a small bed amplitude, one can without loss of generality consider 
sinusoidal perturbations of the form e lkx / H o ; where k = kHo is the wavenumber rescaled 
by the flow height. We note qi, Z%, Ri and Hi the respective disturbances of q, Z, R 
and H. At the linear order, we get: 

Z\ + i?i = H\ 
3i = 

-ikFr 2 R l = tsmORi - ikHi + 2t&n6R l + -tan0P#i (2.17) 

Combining the first and last equations, we obtain the deformation of the free surface as 
a function of that of the bed. 

= 3tanfl + zfcFr 2 

1 (3 + ifc 2 )tan6» + ifc(Fr 2 - 1) K ' ' 

In figure [21 we plot the amplitude and the phase of the free surface profile with respect 
to the sinusoidal bed deformations as functions of the Froude number. Due to viscous 
damping, the bed and the free surface are out of phase essentially close to a unit Froude 
number, and the amplitude of the perturbation is strongly reduced. In contrary, for large 
Froude numbers, the disturbances of the free surface profile exactly match that of the 
bed profile in phase as well as in amplitude. 
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Figure 3. Temporal growth rate a (a) and phase velocity c (b) as functions of the wavenumber 
k for tan# = 0.05 and L sa , t /Ho = 0.06 at three different Froude numbers: Fr = 1.3, Fr — Fr c 
and Fr = 1.4 (thin solid, bold dashed and bold solid lines respectively). 



2.5. Application: anti- dunes in a viscous stream 

We can use the previous calculation to predict the formation of sand dunes in a viscous 
stream. The following calculation is only an illustration of the integral method. We thus 
restrict ourselves to exhibit the interest of the new equations in comparison to the lubrica- 
tion approximation. We refer the reader interested in dunes generated by a steady viscous 
nowtolYaliriT977llYalin 19851 IKuru et al 19951 [Coleman fc Melville 1996llRaudkivi 19971 
Coleman fc Eling 2000(1 Valance fc Langlois 2005l|Kouakou fc Lagree 20051 EE arru &: Hinch 2006 , 



Charru 20061 The linear stability analysis of a flat sand bed can be described in the frame- 



work developed bv lAndreotti et al. 20021 lElbelrhiti et al. 20051 [Claudin fc Andreotti 20061 
for dunes forming in a turbulent boundary layer. The evolution of the bed Z is governed 
by the sand mass conservation: 

d t Z + d x Q = 0, (2.19) 

where Q is the volumic sediment flux, i.e. the volume of the sand grains per unit time 
which cross an infinite normal surface of unit transverse width. The sand flux Q itself 
depends on the basal shear stress r, with a space lag L sat ([Andreotti et al. 2002[> 

isat d x Q = Q sat - Q, (2.20) 

where the saturated sand flux Q sat is a function of r only. For the following analysis, we 
do not need to know the details of the transport law Qsat (t) but only the value of this 
saturated flux over a flat bed, Qq, and that of the coefficient: 

[i= \<m^ (221) 

Qo dr 

The unstable modes are now of the form e (< T -^) t + ikx / H t\ a j s the temporal growth 
rate of the dunes and c = u/k is their phase velocity. Equation (|2 . 19[) then reads 
(a — ito)Zi + ik/HoQi — 0. Similarly, equation (|2.20p simplifies into ikj 'Ho-^satQi = 
PQoTi — Qi, where subscripts 1 denote first order corrections. Combining these two 
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0.1 1 10 100 1000 k 

Figure 4. Diagram of stability for such instability, where Fr is plotted as a function of k for 
the following parameter values: tan 8 — 0.05 and L sa .t/Ho = 0.06. Above the line, the growth 
rate is positive and the perturbation is unstable. The dashed line represents the most unstable 
mode. The thin line corresponds to the very same equations but without the terms coming from 
(|2.12[) . In this case, there is no upper cut-off for k and the most unstable mode is at k — > oo. 



equations, we get 

ikQoP ti 



(2.22) 



(H a + ikL sat ) Z x 

which relates a to the ratio t\/Z\, governed by the hydrodynamics. We can see as from 
now that the sand bed instability can be related to the phase of the shear stress with 
respect to the bed (jElbelrhiti et al. 2005| IClaudin fc Andreotti 2006]) . 

To determine Tx/Z\, we assume that the time scales over which the bed and the flow 
evolve are well separated. For the hydrodynamics, the sand bed can thus be considered as 
static, as in the previous section. One can compute the basal shear stress from equation 
(|2.7|) with t = v d z u\ z=Q = 2>vq/R 2 . Its first order perturbation t\ then reads: 

„ . . „ tan#fc 2 — 2ik „ ,„„„.. 

n - -2 9B *0 Rl =9^e i3+1 _ h2)tan9 + fk{pr2 _ i) Zl . (2.23) 

Inserting this hydrodynamical relation into the expression of the growth rate (I2.22[) . we 
finally obtain: 



Ho , : \ _ - 
— — —{a -iuj) = o-iuj 
) pg sm9 



(2 + it&n9k)k 2 



(1 + ikL sat /H ) [(3 + \k 2 ) tan6» + ik(Fr 2 - 1)] ' 



(2.24) 



where a and uj are the dimcnsionless growth rate and frequency as defined above. We 
also write c = ui/k. 

Both a and c are represented as a function of k in figure [3] for different values of the 
Froude number. For Fr < Fr c , the perturbation remains stable for all wavelengths (a 
is negative). Above Fr c , the perturbation is unstable for a range of k. These lower and 
upper cutoffs are represented in the stability diagram of figure H] Note that the value 
of Fr c depends on the other parameters tan# and L sat /H . Unlike the growth rate, 
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the velocity of the waves is always negative (opposite to the fluid one) below or above 
the critical Froude number. As these sand dunes propagate up-stream, they correspond 
in fact to 'antidunes'. Performing the very same calculation without the terms coming 
from (|2.12p , one get a similar scenario as that displayed in figure [3] with the important 
difference that when perturbations are unstable (above Fr c = 1 independently of the 
other parameters) , the corresponding range of wavenumbers extends up to infinitely large 
k (no upper cut-off). As a conclusion, the longitudinal diffusion of momentum plays a 
crucial role in the control of this instability. Our equations are thus able to recover the 
correct dynamical mechanisms that were present in the original 3D equations but were 
lost in the shallow water approximation - and a fortiori in the lubrication approximation. 



3. Turbulent flow: the Prandtl mixing length approach 

In the previous section, we have shown how it is possible to extend Saint- Venant equa- 
tions beyond the classical shallow water approximation. The method was illustrated on 
the extensively studied viscous case, which extends the lubrication equations to slightly 
incrtial flows. We now apply the same method to fully turbulent flows, for which equation 
(|2.6p strongly under-estimates the rate of energy dissipation. 



3.1. Reynolds stress and Prandtl mixing length closure 

At large Reynolds numbers fluctuations can no longer be neglected. Writing the Navier- 
Stokes equation for the mean flow (intended as an ensemble average over realizations), a 
turbulent stress Tjj = — (u '^u'j) adds to the viscous one, called the Reynolds stress tensor 
(jLandau fc Lifschitz 1959). For large values of Re — UH/u, the viscosity term may be 
neglected and the Saint- Venant equation becomes: 

[dtu + ud x u + wd z u] dz = / [gsinO — d x p + d x r xx + d z r xz ] dz. (3.1) 



In the case of a homogeneous turbulent boundary layer, Tij can be related to the velocity 
gradient using dimensional analysis. The only length-scale is given by the distance to the 
boundary z and the only time scale by the inverse of the shear rate The shear stress 
is then of the form £ 2 |7|7, where I — k(z + zq) is the Prandtl mixing length, k ~ 0.4 
is the von Karman constant, determined experimentally, and Zq is the hydrodynamic 
roughness. For a homogeneous shear stress t xz , one obtains a logarithmic velocity profile 



(3.2) 

In the general case, there exists no univoque relation between the Reynolds stress Tij 
and the velocity field. However, situations close to the homogeneous turbulent bound- 
ary layer can be reasonably described by a phenomenological closure. For the sake of 
simplicity, we have chosen here to introduce a turbulent viscosity v t , derived from the 
Prandtl mixing length approach. The Reynolds stress tensor is written as a function of 
the velocity gradient: — £ 2 \djUi\djUi. In the case of a river, it is commonly assumed 
that the pressure field is hydrostatic and that the velocity profile, as in the case of the 
turbulent boundary layer, is logarithmic. We thus write 

--^Imfi + ^V (3.3) 



K R \ Zq 

p = Patm+ g{H - z)cos6, (3.4) 
where q is the instantaneous local flow rate defined in (|2.5p . The vertical integration of 
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this profile u(z) gives for self-consistency 



kR 



(R + go) In (l 



- R 



(3.5) 



This so-called Chezy coefficient is an important quantity as it gives the ratio of the shear 
stress on the bed to the square of the mean velocity: t xz (z = 0) = C z (q/R) 2 (see below). 

In the steady case, due to the linearity of the shear stress with z, we get the following 
expression for the Prandtl mixing length: 



£ = k(z — Z + Zq) 



H - z 



R 



(3.6) 



We recover, with this phenomenological expression, that the size of the mixing eddies is 
limited by the distance to the rigid boundary as well as by that to the free surface. By 
analogy with the previous viscous equations, we introduce a turbulent viscosity v t in the 
definition of the shear stress tensor Ty = v t djUi. In contrast to the viscous case where 
the viscosity is a constant, v t depends on x, z and t: 



v t = P\d z u\ = Ky/Cz q 



(z 



zo)(H 



R 2 



(3.7) 



3.2. Integral method 

We consider a 2D turbulent fluid flow down an incline. Keeping the relations (|3.3[) and 
(|3.4p . we can express the components of the stress tensor as follows: 



Or 



c z 



c 

2 H 



Zq) In 1 



Z 



c 2 



H 



R 



R 



(3.8) 



Summing up all contributions, the turbulent Saint- Venant equation (|2.14[) reads 
R + Z = H, 
d t R + d x q = 0, 



<>i'l + '>, ( " 77 j = gR cos 9 (tan 9 - d x H) - C z -^ - i 



R 



R 2 (5/2-3X)d x 



q 2 
Cz R 2 



18C z j^[Rd xx Z + 2{d x Zf + d x R8 x Z] 



R 2 



2R[9d x Z +{l~3X)d x R]d x C. 



(3.9) 



with A = In 1 



In the long wavelength approximation, one can neglect in the above equation the 
terms corresponding to longitudinal 'diffusion' of momentum. Furthermore, under the 
approximation z -C R, the value ofa~ 1 + 1/(A — l) 2 is close to unity. Therefore, one 
recovers the standard shallow water equation: 



d t q + d x 



gR cos 9 (tan 9 - d x H) - C z 



R 2 



(3.10) 



Among the additional terms, one recognizes a term corresponding to a diffusion of mo- 
mentum along x. As in the viscous case, the other dissipative terms originate from the 
geometry of the flow. 
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3.3. Example: transverse velocity profile in a river 

To illustrate this method in the turbulent case, we study the transverse mean velocity 
profile in a river homogeneous along the streamwise direction. This is a standard hydro- 
logical geometry, usually treated within the shallow water approximation (Julicn 2002). 
We show here that the transverse momentum fluxes plays a significant role. In other 
words, the predicted mean flow is sensitive to the way friction on the bed and the banks 
is taken into account (jFourriere et al. 2 007). As previously, we call x the direction of 
the flow, tan# the slope of the river and z the axis normal to the bed. y designates the 
direction perpendicular to the flow. For this example, we take a simple polynomial form 
for the transverse depth profile R(y), shown in figure O^a). Such a flat and elongated 
shape is typically that of channels or river beds. Note that the above equations have 
been derived for a flow over a bottom varying along the streamwise direction x in order 
to allow for a straightforward comparison between the three cases (viscous, turbulent 
and granular flow rules). For the determination of the velocity profile in a river, the bed 
varies along the direction transverse to the flow, y. It thus requires a specific derivation 
of the integral method. As it strictly follows the same procedure, we shall mainly explicit 
the main steps and emphasise the differences. 

At first sight, such a channel geometry is close to the homogeneous turbulent boundary 
layer. This suggests that a Prandtl-like closure, although phenomenological, may be 
successful. However, four different distances come into play to construct a mixing length 
I: the distances to the bed and to the free surface as in the previous sub-sections (i.e. in 
the z direction), but also the distances to the left and right banks of the river (in the y 
direction). For the sake of simplicity, we decouple the z and y contributions by writing 



l 2 z \d z u\d z u with £ z = k(z — Z + zq) 
l 2 \dyu\d y u with l v = n'(Sy + y ) 



H 



R 



and 



1 

6y 



W 



I) 



(3.11) 
(3.12) 



In the above expression, yo is the roughness of the river banks - we shall below simply take 
yo = zq as well as k' — n. The choice of a so-called 'harmonic mean' for the computation 
of 5y is simply a trick to get, in a smooth way, the smaller (and thus limiting) value 
among the distances to the banks y and W — y. In this case, the water free surface is 
flat and is taken as the reference, so that H = 0, i.e. Z(y) = —R{y). Keeping again 
the relations (|3 . 3[) and (|3.4|) for the velocity and pressure profiles, one can then explicitly 
compute all stress components 



= c z \u\u 



H 



R 



($y + yo) 2 



2 In 1 



In" 1 



Z 



dyZ 



ZO 



d y (VC z U) d y (Vc z U) 
dy(^C~ z U) ^fC~ z U + C z 



dyZ 



z - Z + z 



(3.13) 



(3.14) 



where we have introduced the notation U = q/R, so that the Saint- Venant equation 
reads: 



8 t q = gRsm6- C Z \U\U + 8 y {(Sy + y ) 2 



CR 



d y (VC z U) d y {^C z U) 



+ A 



dy(VC z U) 



R 



dyR^C z U+ , 

z (R + z a ) 



\dyR\dyRC Z \U\U 
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-1000 - 1 

Figure 5. Transverse profiles of (a) the river bed, (b) the mean flow velocity and (c) the two 
components t xz (long dashed line) and r xy (solid line) of the shear stress on the bed. Lengths 
are in units of the bed roughness zq — y . The velocity is rescaled by \fgz~o sin 9/k, and the 

stresses by gzo sin 0/k. For this example, we choose R(y) = D (4 y/W) 10 (1 — y/W) 10 with 
channel deepest depth D = 1000zo and channel width W = WD. In panel (b), the solid line 
corresponds to the velocity field computed taking into account the friction due to the banks. 
It is qualitatively different and quantitatively lower than that including the bed friction only 
(dashed line). 



Sy + Vo 



C z \d v R\ 



\U\U, 



(3.15) 



with C = (1 + z /R){\ 2 - 2A) + 2. 

This rather complicated equation can be greatly simplified under the assumptions that 
\d y R\ <C 1 (weak transverse slope) and z <C R, which gives at leading order 

8 t U = gsinO- C z ^+d y j*: 2 [y + y (l - ^)] ' \8 V U\ d y ll} . (3.16) 
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It is important to note that the same expression is obtained if the logarithmic velocity 
profile is approximated by a plug flow, under the assumption zq -C R only - it is a 
standard starting point in hydraulics. Whether or not the logarithmic profile is a better 
test function than the uniform one yet remains an open question that would require a 
specific study. 

In the above expression, we can directly read that the last term is of diffusive form, 
with a turbulent eddy viscosity written in terms of £ y and the mean velocity U gradient. 
Without its contribution, i.e. without taking into account the friction of the banks, the 
steady solution of this equation is the usual hydrological relation (Julien 2002) 



Uo = f-^- (3.17) 

The result of the integration of equation (|3.16p is displayed in figure E^b) with the solid 
line. The important point is that, although the typical bed profile curvature scales on the 
water depth, that of the velocity is more like the channel width. This is in agreement with 
experimental and field measures (Fo urriere et al. 2007|) . For comparison, the dashed line 
corresponds to the velocity computed with formula (|3.17[) . Its profile is, by contrast, 
as flat as R(y). Finally, the stress profiles are shown in panel (c). As expected, the 
curvature of the curve T xz (y) is similar to that of U(y). More interesting is the profile 
T xy (y): it goes very quickly (almost at the bank) to its maximum (absolute) value and 
then gently decreases down to zero at the center of the channel. This large value of r xy 
is of the order of that of t xz at the center of the channel where the flow velocity is the 
largest. This means that the banks are in fact as stressed as the center of the channel 
bed. This remark may be of importance as far as erosion processes are concerned. 



4. Dense granular flows 

The integral method can be applied to the case of a dense granular material flowing 
down an inclined plane, just like for a fluid. As for the above turbulent flow case, the 
starting equations are identical to (|2.ip but the viscous term isAu is replaced by djT X j. 
Of course, the constitutive relation, i.e. the relation between the stress tensor to the 
velocity field, has to be specified. We use here the recent advances in this field obtained 
bv lGDR MIDI 20041 IdaCruz et al 20051 |Jop et al. 2006| 

4.1. Local rheology 

It has been shown that the rheology of dense granular flows is local in first approximation. 
This means that the stress and the shear rate at a given location in the flow are related 
through a univoque relation. From dimensional analysis the rheological law can be 
written as rjp — p,{I), where the friction coefficient fi is a function of the properly 
rescaled shear rate I =\j\d/^/p (|GDR MIDI 2004P and d is the grain diameter. Following 



Jop et al. 2006[ we split the stress tensor as —pSij +ry, where p is the isotropic pressure, 



and write the friction law for a granular material as 



^ lr 3 (4.1) 



where jij = d{Uj + djUi is the strain tensor, (7! = y ^ijjij is the second invariant of 
7ij- 
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4.2. Steady uniform solutions 

For the sake of simplicity, we shall restrict to a 2D geometry and look for the homogeneous 
and steady solutions. The momentum balance reads 

= gsin6 + d z T xz . (4.2) 

As the shear stress vanishes at the free surface of the dense granular flow, we get 

t xz = /J.(I)p = g(H Q - z) sin 0, (4.3) 

where Hq is the depth of the flowing layer. Similarly, the isotropic pressure increases 
linearly with depth: p = g(Ho — z) cos0. The ratio of the two gives a constant friction 
coefficient: /i = tan#. As all other components of the strain tensor vanish, the shear rate 
j xz is then selected by the relationship: 

/o = %^=/^ 1 (tan0), (4.4) 

where denotes the inverse function of \x. The integration along the z-axis of j xz 
leads to a Bagnold-like velocity profile 

H*' 2 - (H„ - zfl-' 



Note that this profile is in good agreement with recent two-dimensional numerical simu- 
lations of sheared granular layers (|da Cruz et al. 2 005). 

4.3. Integral method 

At the sight of the above homogeneous steady solutions, we introduce the following test 
function for u 

and keep a hydrostatic pressure profile: 

p = Patm + g(H - z) cos 0. (4.7) 

B, H and q denote the same quantities as in the previous sections. As for ry, we first 
need to compute jij. With the above expression of u, one could in principle compute w 
from the continuity equation. Its expression is rather heavy, but under the assumption 
u w we get to the first order in w/u: j xx — —j Z z = 2d x u and j xz — d z u, so that 
the second invariant of the strain tensor reduces to | --y j = d z u. Finally, the stress tensor 
components read: 



fi(I) g cos 



— y/R(H-z)(Rd x q - qd x B) 



- 2(H - z)d x H + - z f [^ q d x B-Bd x q 



(4.8) 

n(I)g(H - z) cos 0. (4.9) 



The effective friction coefficient /i(/) can be determined experimentally using homoge- 
neous steady flows (|GDR MIDI 20041 IdaCruz et al. 20"05l|Jop et al. 2006[). While sev- 



eral forms have been proposed for fx(I) ( |Pouliquen 1999a Ida Cruz et al. 20051 Jop et al. 2006 
lAndre otti 2007J), the linear expansion around its static value tan#, 

H(I) =taa0 + m(I - I ), (4.10) 
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is sufficient for our practical purpose. For the glass beads used by |Pouliquen 19 99a 
(d = 500 /im), the coefficient m = dfi/dl is around 0.5. Within this linear approximation, 
the Saint- Venant equation gives: 



d t q + d x [a 



R 



-gR cos 6[m{I - I ) + d x H] 



H 



(4.11) 



where a = 5/4 due to the Bagnold-like profile for u. Making the last term explicit, we 
can write the full set of equations for the dense granular case as 



Z + R 
d t R + d x q 



H 




d t q + d x [a 



R 



gRcos9 I m(I — To) + d x H + tan( 
-R 2 d xx q (^(d x Rf + 28 x Hd x Z 



R[ 7 -d xx R + d xx Z 



^ 1 R 2 {d x q) 2 -^-Rd x qd x R 
iq A 3q 



(4.12) 



4.4. Example: the Kapitza instability 

Several experimental works report that surge waves deform the free-surface of granular 
flows down a rough inc lined plane (|Savage 19791 IDaerr 200T| | Fort err e fc Pouliquen 2003j 
Malloggi et al. 20061 . |Forterre fc Pouliquen 2003| has shown that this instability is of the 
same nature as that observed in classical fluids (hence the name 'Kapitza instability') but 
with specificities due to the granular rheology. We show here that our equations allow 
to recover the correct characteristics of this instability, lost with the long wavelength 
expansion derived by |Forterre k, Pouliquen 2 003 



As the bed is now flat, the first equation of (|4.12|) simply reduces to R = H . In the 
limit of small perturbations, we can look for a solution of the form: 



H = H ( 1 + H l e t{hx/Ha ' atqo/H » ) 



and 



go U + qie 



i(kx/H -Otqo/H%)) 



(4.13) 



where k = kHo denotes the rescaled wavenumber and lu = ojHq / qo the dimensionless 
angular frequency. To the first order, the system of equation (|4.12jl becomes: 



Fr l 



— iuHi + ikqi 
2ik I qi — ikHi 



= 0, 

= y-t&nOk 2 



-mlo — ik I Hi 



mlo 



where the Froude number of the flow is defined by: 



4 - 9 
-tan 9k 2 
9 



Fr z = a 



gHfi cos 6 



k and w are finally related by the implicit dispersion relation: 
1 9 z- mI .5 T . .1 _-9 tan 9 



(4w - 7k)k 2 



(4.14) 



(4.15) 



(4.16) 



By contrast to the dispersion relation of ( jForterre fc Pouliquen 2003| , the above one 
now includes new terms - those gathered on the right hand side. This instability being 
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u i 1 1 1 1 

1 CO 2 

Figure 6. Theoretical (lines) and experimental (•) dispersion relation for 9 = 29° and 
Fr = 1.14. Spatial growth rate s (a) and phase velocity c (b) as functions of the pulsation 
Q. Inset (c) shows c for larger values of Q. Note the fall-off of the curve around lj ~ 20. The 
dashed lines represent theoretical curves of |Forterre fc Pouliquen 2003 1 with a = | — In their 
original paper they rather chose a = 1. Lengths are in units of Ho and times in units of qo/Hq. 

convective, a real w is imposed and the instability develops spatially downflow. We 
then split the wave number into two terms k = Cu/c — is, where c and s represent the 
(dimensionless) phase velocity and the spatial linear growth rate. Finally the functions 
s(u>) and c(u>) are computed from (|4.16[) . 

Equation (|4.16[) is of the third order and among the three solutions, two correspond 
to relaxation modes: the product sc is negative for all ui and all Froude numbers. The 
analysis of the third one shows that the flow is unstable when the Froude number exceeds 
the critical value Fr c — 1 (independently of 0). For Fr > 1, small angular frequencies 
are unstable up to a cut-off value Co c - This cut-off frequency vanishes at Fr — Fr c , which 
means that it is a zero mode instability. The re-stabilization of the flow for frequencies 
larger than uj c is missed without the new terms in the dispersion relation (|4.16[) . 

As for this problem, experimental data ( |Forterre fc Pouliquen 2003] ) as well as the 
results of a numerical integration of the full equations (jForterre 2006|) are available, it is 
interesting to discuss some quantitative aspects of the integral method predictions. The 
results of our calculations as well as the data are displayed in figure [H In panel (a) , 
one can see that, with Fr = 1.14 corresponding to the experiments, the estimation of 
cJ c is very good, however, the magnitude of the maximum growth rate s max is under- 
estimated by a factor of 3. A similar quantitative difference is found for the value of 
the critical Froude number: its experimental determination is Fr c ~ 0.6, i.e. almost 
twice smaller than our prediction. For comparison, the computation of s max and Fr c 
in IForterre 20061 are respectively a factor of 2 below and 10% above the corresponding 
experimental values. In panel (b) is displayed the velocity of the waves. The agreement 
between theory and experiment is fair. Interestingly, our calculation predicts a fall-off 
of the curve around Q ~ 20 and a cross-over to a second plateau, see panel (c). It is 
worth noting that this change of behaviour at large frequency is also predicted in the 
liquid case. This unexpected behaviour of the curve c{ui) could not be checked with the 
calculations of Forterre 2006] as, in accordance with the experimental data, the numerics 
was run up to ui — 2 only. 
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5. Summary and perspectives 

In this paper, we have shown that, for the description of the dynamical properties 
of flows, one can go beyond the lubrication or shallow water approximations following 
a well-defined procedure that can be adapted to any constitutive relation: the integral 
method. The equations derived within this framework are simpler to integrate than 
the full 3D equations (e.g. the Navier-Stokes equations) but possess the same physical 
ingredients. The technique consists in choosing 'test functions' for the velocity profiles, 
whose functional form is inferred from the exact homogeneous and steady solutions of 
the full equations, and in performing a depth averaging procedure with these profiles. 

We have successfully applied this method to situations as diverse as viscous, turbulent 
and granular flows down an incline. In these three cases, we have investigated a specific 
example, respectively, formation of anti-dunes, velocity profile in a river cross-section and 
development of the Kapitza instability. We have emphasised the benefits of the method 
compared to the usual approximations. In particular, we have shown that important 
physical mechanisms are qualitatively recovered: the presence of wavelength or frequency 
cut-offs in the dispersion relations of anti-dune or Kapitza instabilities, the significant 
role of the friction at the river banks on the velocity profile curvature. 

A strong assumption in the derivation proposed concerns the exchanges of vertical mo- 
mentum. In the three cases, we have indeed assumed that the pressure was dominated 
by gravity effects. It has for instance been shown by |Savage fc Flutter 1991] that the to- 
pography induces corrections to the hydrostatics pressure profile. It has also been shown 
by |Mangeney et al. 2005 that, in the first stages of the vertical collapse of a granular 
column, gravity is not balanced by the pressure but by the vertical acceleration. The de- 
scription of such a situation is beyond the scope of the equations derived here. However, 
the same strategy could be applied, using the depth-averaged vertical momentum equa- 
tion to determine the pressure, whose profile would be prescribed by a test-function. By 
analogy to the choice made for the velocity, a linear pressure profile is probably the best 
choice as one recovers the above equations when vertical inertia and vertical diffusion of 
momentum are negligible. A dedicated work is needed to test this idea further on. Our 
aim is now to use the integral method as a tool to investigate physical problems such 
as the formation of current ripples versus dunes in viscous or turbulent flumes, the river 
width selection or the description of fingering effects in granular avalanches. 
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